# -*-   
# Author: qinguoming
# Date: 2025-02-11 
# Version: 1.0 
# Description: 
#      1.  这个程序是一个模板,用于对2D线性桁架单元组成的结构的静态分析,输出节点位移,输出单元应力,单元应变,节点支反力等结果.
#      2.  2D线性桁架节点自由度为: ux,uy      
#      3.  节点和网格的编号必须从1开始,且连续.
#   脚本执行: \EF2D> python -m Scripts.Script_for_2d_Truss_StaticAnalysis
# -------------导入模块------------------
""" 
    from EF2D.Elements import T2D2
    from EF2D.node import Node
    from EF2D.Loads import DloadOnTruss, ConcentratedLoad,DisplacementBc
    from EF2D.Solver import *
    from EF2D.fea import FEA
    from EF2D.mesh import MeshBlock
    from EF2D.tools import *
    from EF2D.InpReader import read_abaqus_inp 
"""
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
from typing import List,Dict
from icecream import ic
import sys
import os
from EF2D import *
# -------------定义变量------------------
job_name="job-1"
# 网格文件
mesh_file="L:\\EF2D\\tests\\test_Script_2d_truss_static.inp"
# 材料属性(SI-m单位制)
rho=7300
E=2.1e11
A=1e-4
# 约束:支持位移边界条件
# 格式:name,node_label_list,dof,value
# 例如:
def_bc=[("fix_ux","fix_x",1,0.0),
        ("fix_uy","fix_y",2,0.0)]
# 定义载荷(载荷类型,载荷参数)
# dload:代表DloadOnTruss; cload:代表ConcentratedLoad
# def_load=[("dload","load_1000N",[1,2,3,4],1000.0),
#             ("cload","load_2000N","nset_name",2,2000.0)]
def_load=[("cload","load_-1000N","load_y_n1000N",2,-1000.0),]

# ////////求解设置//////////////////////
bcMethod="default"  # 边界条件处理方法,默认使用充大数法; 其他值为1,2,3,4分别是Solver.py的方法
export_result=True  # 是否输出结果
export_file="result.mat"  # 输出文件名
# -------------导入网格------------------
node_data, element_data, nset_data, elset_data ,_= read_abaqus_inp(mesh_file)
print("节点数:",node_data.shape[0])
print("单元数:",element_data.shape[0])
print("节点集:",nset_data.keys())
print("单元集:",elset_data.keys())
print("边界条件:",[b[0] for b in def_bc])
print("载荷有:",[l[1] for l in def_load])
print("边界条件处理方法:",bcMethod)

# -------------建立网格------------------
node_num=node_data.shape[0]
element_num=element_data.shape[0]
# 定义节点对象
node_dict:Dict[int,Node]={}
for i in range(node_num):
    if node_data[i,0] not in node_dict.keys():
        node_dict[node_data[i,0]]=Node(label=node_data[i,0],x=node_data[i,1],y=node_data[i,2],z=node_data[i,3])
# 定义单元对象
element_list:List[T2D2]=[]
for i in range(element_num):
    element_list.append(T2D2(label=element_data[i,0],
                                        node1=node_dict[element_data[i,1]],
                                        node2=node_dict[element_data[i,2]],
                                        E=E,A=A,density=rho))
# 定义MeshBlock对象
block_list=[MeshBlock(name="block1",elements=element_list),]
# 定义载荷和边界
bc_list=[]
for b in def_bc:
    if isinstance(b[1],str) and b[1].lower() in nset_data.keys():
        bc_list.append(DisplacementBc(name=b[0],nodeLabels=nset_data[b[1].lower()],dofth=b[2],value=b[3]))
    else:
        bc_list.append(DisplacementBc(*b))
load_list=[]
for l in def_load:
    if l[0]=="dload":
        if isinstance(l[2],str) and l[2].lower() in elset_data.keys():
            load_list.append(DloadOnTruss(name=l[1],elemLabels=elset_data[l[2].lower()],value=l[3]))
        else:
            load_list.append(DloadOnTruss(*l[1:]))
    elif l[0]=="cload":
        if isinstance(l[2],str) and l[2].lower() in nset_data.keys():
            load_list.append(ConcentratedLoad(name=l[1],nodeLabels=nset_data[l[2].lower()],dofth=l[3],value=l[4]))
        else:
            load_list.append(ConcentratedLoad(*l[1:]))

# -------------建立FEA模型-------------------------
Struc=FEA(meshblocks=block_list,bcs=bc_list,loads=load_list)

# -------------运行分析------------------
if bcMethod=="default":
    ug,fg=Struc.Solve(StepType="Static")
else:
    Kg=Struc.Kg()
    bc_dict=Struc.ProcessBcs()[-1]
    load_dict=Struc.PreProcessFg()[-1]
    if bcMethod==1:
        ug,fg=Method1(Kg,bc_dict,load_dict)
    if bcMethod==2:
        ug,fg=Method2(Kg,bc_dict,load_dict)
    if bcMethod==3:
        ug,fg=Method3(Kg,bc_dict,load_dict)
    if bcMethod==4:
        ug,fg=Method4(Kg,bc_dict,load_dict)
Struc.Ug=ug
Struc.F_react=fg

# -------------输出结果--------------------
if export_result:
    el_stress=np.zeros((element_num,1))
    el_strain=np.zeros((element_num,1))
    # 计算应力、应变、节点位移、节点支反力
    for i,e in enumerate(element_list):
        el_stress[i,0]=e.cal_stress(ug[e.ElemDofIndexs,:])
        el_strain[i,0]=e.cal_strain(ug[e.ElemDofIndexs,:])
    # 输出结果到mat文件,可以和matlab对接
    mat_data={
        "Node_data":node_data,
        "connecticity_data":element_data,
        "vtk_type":element_list[0].VTK_TYPE,
        "node_num":int(node_num),
        "element_num":int(element_num),
        "node_dof_num":2,
        "node_ux":ug[::2,0],
        "node_uy":ug[1::2,0],
        "node_u_mag":np.sqrt(ug[::2,0]**2+ug[1::2,0]**2),
        "node_Freact_x":fg[::2,0],
        "node_Freact_y":fg[1::2,0],
        "node_Freact_mag":np.sqrt(fg[::2,0]**2+fg[1::2,0]**2),
        "element_stress":el_stress,
        "element_strain":el_strain,
        }
    sio.savemat(export_file,mat_data,appendmat=False,oned_as='row')
    print("结果已输出到文件:",os.getcwd()+os.sep+export_file)
